### This script will compare empirical measurements of microclimate to three climate predictions
### The three predictions are BIOCLIM, AWAPCLIM, and microCLIM
### Should also plot daily microclimate for daily microCLIM predictions

### Load libraries

library('SDMTools')

### Set working directories

out.dir = '/home1/99/jc152199/ChapterOne/'

### Read in Raw Data file with preds

raw.data = read.csv('/home1/99/jc152199/brt/FINALOPTIMALMODELS/Model_Data_Plus_Preds.csv',header=T)


######################### Jump down to next multi hash line

outdata = NULL

for (mod in c('max','min')) # Loop through both models, determining the 95% CI's for slope and intercept of both lm's

	{
	
	t.data = data.frame(site=raw.data$site,raw.data[,grep(paste(mod),names(raw.data))]) # Subset data by surface type
	
	AWAP_lm = lm(t.data[,2]~t.data[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	mc_lm = lm(t.data[,2]~t.data[,4]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	
	t.AWAP.sigdata = data.frame(pred='AWAP',mod=mod, min.slope = summary(AWAP_lm)$coefficients[2,1]-(1.96*summary(AWAP_lm)$coefficients[2,2]), slope.pred = summary(AWAP_lm)$coefficients[2,1], max.slope = summary(AWAP_lm)$coefficients[2,1]+(1.96*summary(AWAP_lm)$coefficients[2,2]), min.int = summary(AWAP_lm)$coefficients[1,1]-(1.96*summary(AWAP_lm)$coefficients[1,2]), int.pred = summary(AWAP_lm)$coefficients[1,1], max.int = summary(AWAP_lm)$coefficients[1,1]+(1.96*summary(AWAP_lm)$coefficients[1,2]))
	t.mc.sigdata = data.frame(pred='mc',mod=mod, min.slope = summary(mc_lm)$coefficients[2,1]-(1.96*summary(mc_lm)$coefficients[2,2]), slope.pred = summary(mc_lm)$coefficients[2,1], max.slope = summary(mc_lm)$coefficients[2,1]+(1.96*summary(mc_lm)$coefficients[2,2]), min.int = summary(mc_lm)$coefficients[1,1]-(1.96*summary(mc_lm)$coefficients[1,2]), int.pred = summary(mc_lm)$coefficients[1,1], max.int = summary(mc_lm)$coefficients[1,1]+(1.96*summary(mc_lm)$coefficients[1,2]))

	outdata = rbind(outdata,t.AWAP.sigdata,t.mc.sigdata)
	
	}
	
# Close

# Determine signicance of slopes (must fall between the 95% CI's to be a non-significant difference

outdata$slope.sig = NA
outdata$slope.sig[which(outdata$min.slope>1)]='sig'
outdata$slope.sig[which(outdata$max.slope<1)]='sig'
outdata$slope.sig[which(outdata$min.slope<1 & outdata$max.slope>1)]='nonsig'

# Determine signicance of intercepts (must fall between the 95% CI's to be a non-significant difference

outdata$int.sig = NA
outdata$int.sig[which(outdata$min.int>0)] = 'sig'
outdata$int.sig[which(outdata$max.int<0)] = 'sig'
outdata$int.sig[which(outdata$min.int<0 & outdata$max.int>0)] = 'nonsig'


# Write out summary data as Table 1

write.csv(outdata, file=paste(out.dir,'Tables/Table1.csv',sep=''), row.names=F)

############################################################################### Start again here





for (mod in c('max','min')) # Loop through both models, overlaying scatterplots of AWAP vs Emp and microCLIM vs Emp

	{
	
	t.data = data.frame(site=raw.data$site,raw.data[,grep(paste(mod),names(raw.data))]) # Subset data by surface type
	
	AWAP_lm = lm(t.data[,3]~t.data[,2]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	mc_lm = lm(t.data[,4]~t.data[,2]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	
	lims = range(c(t.data[,2],t.data[,3],t.data[,4]),na.rm=T)
	
	png(paste(out.dir,'Figures/Figure1',mod,'.png',sep=''), height=18, width=18,units='cm', res=1000) # Open .png device driver	
	
	plot(t.data[,3],t.data[,2],xlab='Empirical',ylab='Predicted',main = paste('Daily T',mod,sep=''),cex=.7, type='n', xlim=lims, ylim=lims) # Configure plot space
	
	#abline(b=1,a=0,col='black') # 1:1 line
	
	#polygon(AWAP.polyout[,1],AWAP.polyout[,2], col='#FF000025', lty=1, lwd=1, border=NA) # Plot a polygon representing the 95% CI's for the AWAP linear model
	
	points(t.data[,2],t.data[,3],col='#FF000050',cex=.8,pch=16) # Points for AWAP and Empirical
	
	points(t.data[,2],t.data[,4],col='#0000FF50', cex=.8,pch=16)
	
	abline(b=summary(AWAP_lm)[[4]][2],a=summary(AWAP_lm)[[4]][1],col='red',cex=2, lwd=2,lty=1) # Plot a line representing the relationship between AWAP and Empirical data
	
	#polygon(mc.polyout[,1],mc.polyout[,2], col='#0000FF25', lty=1, lwd=1, border=NA) # Plot a polygon representing the 95% CI's for the microCLIM linear model
	
	abline(b=summary(mc_lm)[[4]][2],a=summary(mc_lm)[[4]][1],col='blue',cex=2,lwd=2,lty=1) # Plot a line representing the relationship between microCLIM and Empirical data
		
	legend('topleft',legend = c('BRT Model',paste('Adj. r^2 ',substr(summary.lm(mc_lm)[9],1,4),sep=''),paste('Slope ',round(mc_lm$coefficients[2],2),sep=''),paste('Intercept ',round(mc_lm$coefficients[1],2),sep='')),text.col=c('black','blue','blue','blue'), bty='n', cex=.8)
	legend('bottomright',legend = c('AWAP Model',paste('Adj. r^2 ',substr(summary.lm(AWAP_lm)[9],1,4),sep=''),paste('Slope ',round(AWAP_lm$coefficients[2],2),sep=''),paste('Intercept ',round(AWAP_lm$coefficients[1],2),sep='')),text.col=c('black','red','red','red'), bty='n', cex=.8)
    
	
	dev.off()
	
	}
	

	
#### Below this are lines to create a polygon shape to plot

#newdata1 = data.frame(obs=t.final.out[,4]) # Bind x-data for plots into a dataframe
	#newdata2 = data.frame(obs=t.final.out[,3])
	
	#AWAP.rawpoly = predict(AWAP_lm, newdata1, interval="confidence"); AWAP.rawpoly = cbind(newdata1,AWAP.rawpoly) # Calculate 95% CI's for all x values in the AWAP linear model
	#mc.rawpoly = predict(mc_lm, newdata2, interval="confidence") ; mc.rawpoly = cbind(newdata2,mc.rawpoly)# Calculate 95% CI's for all x values in the microCLIM linear model
	
	#AWAP.polybottom = AWAP.rawpoly[c(order(newdata1)),c(1,3)]
	#names(AWAP.polybottom) = c('x','y')
	#AWAP.polytop = AWAP.rawpoly[c(rev(order(newdata1))),c(1,4)]
	#names(AWAP.polytop) = c('x','y')
	#AWAP.polyout = rbind(AWAP.polybottom,AWAP.polytop) # Calculate a polygon that represents the 95% CI's of the AWAP linear model
	
	#mc.polybottom = mc.rawpoly[c(order(newdata2)),c(1,3)]
	#names(mc.polybottom) = c('x','y')
	#mc.polytop = mc.rawpoly[c(rev(order(newdata2))),c(1,4)]
	#names(mc.polytop) = c('x','y')
	#mc.polyout = rbind(mc.polybottom,mc.polytop) # Calculate a polygon that represents the 95% CI's of the microCLIM linear model
	
	#lims = range(c(AWAP.polyout[,1],AWAP.polyout[,2],mc.polyout[,1],mc.polyout[,2])) # Calculate limits for x and y axes
	
	
